*! percentile as estimation procedure
* v.1.2 by Stas Kolenikov, skolenik at gmail dot com

program define epctile, eclass sortpreserve properties(svyb svyj)
  * survey properties of the program?

  version 11
  if !replay() {
     syntax varname [if] [in] [pw fw iw/], Percentiles(numlist) ///
         [ svy over(varlist) Level(cilevel) subpop(string) VALuemask( string ) SPEClabel noSE ]

     * the main study variable
     local y `varlist'

     * sample conditions
     marksample touse
     sort `touse' `over' `y'

     * parse the subpop option
     if "`subpop'"!="" {
        ParseSubPop `subpop', touse(`touse')
        local subpop subpop( `subpop' )
     }

     * parse the weights
     tempvar wgt
     if "`weight'" == "" {
       * no weight specified
       qui gen byte `wgt' = `touse'
     }
     else {
       * something was put down as weights
       qui gen double `wgt' = `exp' * `touse'
       local wexp [pw=`wgt']
     }
     * parse the survey option
     if "`svy'"~="" {
       * are there survey settings to fall back on?
       qui svyset
       if "`r(wvar)'" != "" {
         * if there are any weights found in svy settings, replace the weight var
         qui replace `wgt' = `r(wvar)'
       }
       if "`r(settings)'" != ", clear" {
         * if there are any svy settings at all
         local svy svy , `subpop' :
       }
       else {
         * svy option was specified, but no svy settings found
         di as err "Warning: no svy settings found"
         local wexp [pw=`wgt']
       }
       if "`if'`in'"!="" {
         * if/in or subpop?
         qui count if `touse'
         if r(N) != _N di as err "Warning: if and in conditions should be implemented via svy, subpop option"
       }
     }

     qui sum `wgt' if `touse'
     local sumw = r(sum)

     * check if percentiles are within bounds
     tokenize `percentiles'
     local np : word count `percentiles'
     forvalues j=1/`np' {
        local p`j' : word `j' of `percentiles'
        if `p`j'' < 10 local p`j' 0`p`j''
        if `p`j'' <= 0 & `p`j''>=100 {
           di as err "percentiles must be between 0 and 100"
           exit 198
        }
     }

     qui count if `touse'
     local nobs = r(N)

     tempname b bb V
     * 1. get percentiles via _pctile
     if "`over'" == "" {
        _pctile `y' [pw=`wgt'] if `touse', percentiles( `percentiles' )

        mat `b' = J(1,`np',.)
        forvalues j=1/`np' {
           mat `b'[1,`j'] = r(r`j')
           local cnames `cnames' `=subinstr("p`p`j''",".","d",.)'
        }
     }
     else {
        * handle the OVER option here!
        unab overlist : `over'
        * number of variables
        local novervar : word count `overlist'

        FormXLevels `overlist', valuemask( "`valuemask'" )

        mat `b' = J(1,`np'*`s(ncombos)', .)
        * the parameter estimates should run by percentile, then by `over' variables
        * to match the estimation output of -mean-

        forvalues k=1/`s(ncombos)' {
           _pctile `y' [pw=`wgt'] if `touse' `s(if`k')', percentiles( `percentiles' )
           forvalues j=1/`np' {
              mat `b'[1,(`j'-1)*`s(ncombos)' + `k'] = r(r`j')
           }
        }
        forvalues j=1/`np' {
           forvalues k=1/`s(ncombos)' {
              if "`speclabel'" == "" local cnames `cnames' `=subinstr("p`p`j''",".","d",.)':`s(vmlab`k')'
              else                   local cnames `cnames' `=subinstr("p`p`j''",".","d",.)'_`s(vmlab`k')'
           }
        }
     }

     if "`se'" != "nose" {

     * 2. compute the linearization pieces
     if "`over'" == "" {
       forvalues k=1/`np' {
          tempvar d`k'
          qui gen `d`k'' = (  (`y'<`b'[1,`k']) - 0.01*`p`k'' )
          local dlist `dlist' `d`k''
       }
     }
     else {
       qui forvalues k=1/`np' {
          tempvar d`k'
          gen `d`k'' = .
          forvalues j=1/`s(ncombos)' {
            replace `d`k'' = (  (`y'<`b'[1,(`k'-1)*`s(ncombos)' + `j']) - 0.01*`p`k'' ) if `touse' `s(if`j')'
          }
          local dlist `dlist' `d`k''
       }
     }

     * 3. compute the covariance with -svy: mean-
     if "`over'" != "" local overopt over(`over', nolabel)
     qui `svy' mean `dlist' if `touse' `wexp', `overopt'
     local vcetype `e(vcetype)'
     local vce `e(vce)'
     mat `V' = e(V)

     * sreturn values are lost by now
     if "`over'"!= "" FormXLevels `overlist', valuemask( "`valuemask'" )

     * 4. compute the Woodruff estimates of variance
     tempname scale
     mat `scale' = J(1,`np'*`e(N_over)',.)
     forvalues j=1/`np' {
        if "`over'" != "" {
          forvalues k=1/`s(ncombos)' {
             local pplus  = `p`j'' + 200*_se[ `d`j'':`s(lab`k')' ]
             if `pplus' >= 100 local pplus `p`j''
             * have to use one-sided derivative
             local pminus = `p`j'' - 200*_se[ `d`j'':`s(lab`k')' ]
             if `pminus' <= 0 local pminus `p`j''
             * have to use one-sided derivative
             _pctile `y' [pw=`wgt'] if `touse' `s(if`k')', percentiles( `pminus' `pplus' )
             mat `scale'[1,(`j'-1)*`s(ncombos)'+`k'] = 100*( r(r2) - r(r1) )/( `pplus'-`pminus' )
/*
             mat `scale'[1,(`j'-1)*`s(ncombos)'+`k'] = ( r(r2) - r(r1) )/( 4 * _se[ `d`j'':`s(lab`k')' ] )
             if `pplus' == `p`j'' | `pminus' == `p`j'' {
                mat `scale'[1,(`j'-1)*`s(ncombos)'+`k'] = 2 * `scale'[1,(`j'-1)*`s(ncombos)'+`k']
             }
*/
          }
        }
        else {
          local pplus  = `p`j'' + 200*_se[ `d`j'' ]
          if `pplus' >= 100 local pplus `p`j''
          * have to use one-sided derivative
          local pminus = `p`j'' - 200*_se[ `d`j'' ]
          if `pminus' <= 0 local pminus `p`j''
          * have to use one-sided derivative
          _pctile `y' [pw=`wgt'] if `touse', percentiles( `pminus' `pplus' )
          mat `scale'[1,`j'] = 100*( r(r2) - r(r1) )/( `pplus'-`pminus' )
/*
          mat `scale'[1,`j'] = ( r(r2) - r(r1) )/( 4 * _se[ `d`j'' ] )
          if `pplus' == `p`j'' | `pminus' == `p`j'' mat `scale'[1,`j'] = 2 * `scale'[1,`j']
*/
        }
     }

     * 5. combine into the covariance matrix of the percentiles
     mat `V' = diag( `scale' ) * `V' * diag( `scale' )
     }

     else {
       * no standard errors
       mat `V' = J( colsof(`b'), colsof(`b'), 0 )
     }

     * grand finale: nice labels
     mat rownames `V' = `cnames'
     mat colnames `V' = `cnames'

     mat rownames `b' = `y'
     mat colnames `b' = `cnames'

     mat colnames `scale' = `cnames'
     mat rownames `scale' = `y'

     ereturn post `b' `V', obs(`nobs') depname(`y') esample(`touse')
     ereturn matrix scale  = `scale'
     ereturn scalar N_over = `np'
     ereturn local over      `over'
     ereturn local depvar    `y'
     ereturn local vce       `vce'
     ereturn local vcetype   `vcetype'
     ereturn local cmd       epctile
     * set e(cmd) last
  }
  else { // replay
     if "`e(cmd)'"!="epctile" error 301
     syntax [, Level(cilevel)]
  }
  * output any header above the coefficient table
  di _n as text "Percentile estimation"
  ereturn display, level(`level')
  di _n
end



program define ParseSubPop
   * Jeff Pitblado's suggestion: undocumented _svy_subpop marking tool
   syntax [varlist(max=1 default=none)] [if/] , touse(varname)
   * subpop is [varname] [if]
   * target: zero out observations outside of the subpopulation
   * the input variable is typically `touse'

   confirm variable `touse'

   if "`varlist'`if'" == "" {
      di as err "undefined subpopulation"
      exit 198
   }
   qui {
      if "`varlist'" == "" local varlist 1
      if "`if'" == "" local if 1
      replace `touse' = `touse' & (`varlist') & (`if')
   }
end

exit

History:

v. 1.0 -- March 05, 2010. Basic features: percentiles for i.i.d. data and -svy- via the option
v. 1.1 -- March 06, 2010. Options: over, valuemask, speclabel, subpop,
v. 1.2 -- May 25, 2010.   Fixed an issue with standard errors and -over- option


Simulations:

* i.i.d. data
set seed 10101
pro def epctsim
  drop _all
  set obs 1000
  gen y = uniform()
  epctile y , p(25 50 75)
end

simulate _b _se , reps(1000) saving( epctile_iid, replace) : epctsim
sum

* stratified single stage
clear all
set seed 10101
set obs 500000
gen long id = _n
gen byte h = (_n-1)/10000+1
bysort h (id) : gen m = 2+uniform() if _n == 1
bysort h (id) : gen s = 0.5+0.1*uniform() if _n == 1
bysort h (id) : replace m = m[1]
bysort h (id) : replace s = s[1]
gen y = exp( m + rnormal()*s )
sum y, d
scalar pop05 = r(p5)
foreach k in 10 25 50 75 90 95 {
  scalar pop`k' = r(p`k')
}
sort s id
egen strata = group( h )

save strpop, replace

pro def epctsvysim
   use strpop, clear
   gen r1 = uniform()
   gen r2 = uniform()
   sort strata r1 r2
   bysort strata (r1 r2) : keep if _n < 20+strata
   gen samplw = 10000/(20+strata)
   svyset [pw=samplw], strata(strata)
   epctile y, p(5 10 25 50 75 90 95) svy
end

simulate _b _se , reps(1000) saving( epctile_svy, replace) : epctsvysim
foreach k in 05 10 25 50 75 90 95 {
  gen byte reject_l`k' = ( p`k'_b_cons + invnorm( 0.05 )*p`k'_se_cons > scalar( pop`k') )
  gen byte reject_u`k' = ( p`k'_b_cons + invnorm( 0.95 )*p`k'_se_cons < scalar( pop`k') )
}
sum
foreach x of varlist reject* {
  bitest `x' == 0.05
}

* stratified two-stage
clear all
set seed 10201
set obs 50
gen double m = 2+uniform()
gen double s = 0.5+0.2*uniform()
sort s
gen byte strata = _n
gen int Nh = 100 + rnbinomial( 10, 0.2 )
gen byte nh = 2 + (uniform()<0.4) + (uniform()<0.1)
tab nh
expand Nh
gen double y0 = exp( m + s*rnormal() )
gen int PSU = _n
gen int Mhi = 20 + rnbinomial( 10, 0.1 )
sum Mhi
expand Mhi
gen double y = y0 - ( ln(uniform()) + ln(uniform())  + ln(uniform()) + 3 )
bysort strata PSU : gen byte first = (_n==1)

sum y, d
scalar pop05 = r(p5)
foreach k in 10 25 50 75 90 95 {
  scalar pop`k' = r(p`k')
}

compress

save strpop2, replace

pro def epctsvy2sim
   use strpop2, clear
   gen float r1 = uniform() if first
   gen float r2 = uniform() if first
   bysort strata PSU (first) : replace r1 = r1[_N]
   bysort strata PSU (first) : replace r2 = r2[_N]
   bysort strata (r1 r2 first) : gen iPSU = sum( first )
   bysort strata r1 r2 (first) : replace iPSU = iPSU[_N]
   * sample of PSUs
   keep if iPSU <= nh
   * sample of observations
   gen float r3 = uniform()
   gen float r4 = uniform()
   bysort strata iPSU (r3 r4) : keep if _n <= 10
   gen samplw = Mhi/10 * Nh/nh
   svyset iPSU [pw=samplw], strata(strata)
   epctile y, p(5 10 25 50 75 90 95) svy
end

clear
set mem 80m

simulate _b _se , reps(1000) saving( epctile_svy2, replace) : epctsvy2sim
foreach k in 05 10 25 50 75 90 95 {
  gen byte reject_l`k' = ( p`k'_b_cons + invnorm( 0.05 )*p`k'_se_cons > scalar( pop`k') )
  gen byte reject_u`k' = ( p`k'_b_cons + invnorm( 0.95 )*p`k'_se_cons < scalar( pop`k') )
}
sca li
sum
foreach x of varlist reject* {
  bitest `x' == 0.05
}
